Efficient Isoparametric Integration over Arbitrary, Space-Filling Voronoi Polyhedra 

for Electronic-Structure Calculations 



Aftab Alam, 1 ' 3 S. N. Khan, 2 ^ 3 , Brian G. Wilson, 4 , and D. D. Johnson 1 ' 2 ^ 

1 Department of Materials Science and Engineering, University of Illinois, Urbana, IL 61801, USA 

2 Department of Physics, University of Illinois, Urbana, IL 61801, USA 
3 Division of Materials Science and Engineering, Ames Laboratory, Ames, Iowa 50011, USA and 
4 Lawrence Livermore National Laboratory, Livermore, California 94550, USA 

(Dated: January 25, 2013) 

A numerically efficient, accurate, and easily implemented integration scheme over convex Voronoi 
polyhedra (VP) is presented for use in ab-initio electronic-structure calculations. We combine a 
weighted Voronoi tessellation with isoparametric integration via Gauss-Legendre quadratures to 
provide rapidly convergent VP integrals for a variety of integrands, including those with a Coulomb 
singularity. We showcase the capability of our approach by first applying to an analytic charge- 
density model achieving machine-precision accuracy with expected convergence properties in mil- 
liseconds. For contrast, we compare our results to those using shape-functions and show our ap- 
proach is greater than 10 s faster and 10 7 more accurate. A weighted Voronoi tessellation also allows 
for a physics-based partitioning of space that guarantees convex, space-filling VP while reflecting 
accurate atomic size and site charges, as we show within KKR methods applied to Fe-Pd alloys. 

PACS numbers: 02.60.Jh, 71.15.Dx 



A variety of science and engineering research problems 
require multicentered integrals that cannot be solved an- 
alytically due to the complex domains of integration. 
The total energy and potential in any site-centered, 
electronic-structure calculation involves the evaluation of 
three-dimensional integrals over convex Voronoi polyhe- 
dra (VP) P Although a number of numerical integration 
techniques have been proposedpHsl an efficient, accurate, 
reliable and easily implemented scheme is still lacking. 
Prior methods often rely on detailed analysis of sym- 
metry properties of the integration domains and, hence, 
limit their applicability to arbitrary atomic geometries 
and structures. A major continuing need is an integra- 
tion method over space-filling VP that has a high de- 
gree of accuracy with a minimal computational effort and 
that is sufficiently generic so that it can be used in most 
electronic-structure application codes. 

For example, "exact" linear muffin-tin orbital (EMTO) 
methocP uses an approach from Gonis et al. s ^ to over- 
come VP integration issues for the Poisson potential, 
but it is extremely slowly convergent; various KKR- 
based codes, such as the linear-scaling multiple-scattering 
(LSMS), 9 utilizes shape- functions to perform VP inte- 
grations, which, as we show, is slowly convergent and 
limited in accuracy; the full-potential linear augmented 
wave^ (FLAPW) method avoids VP integrals (via non- 
overlapping muffin-tins and Fourier methods over entire 
unit cell) , but never determining site- VP-specific proper- 
ties and requiring a larger number of spherical harmonic 
basis functions and huge number of plane-waves. 

We present s uch an algorithm by combining a weighted 
VP tessellatiorpJiiJ with isoparametric methods^ to pro- 
vide rapidly convergent integrals for various integrands, 
including Coulomb singularities. For generality, we 
use a Radical Plane Construction^ (RPC) or Power 
DiagraimfS] to guarantee convex, space-filling VP. For 



electronic-structure use, physics-based weights are opti- 
mally chosen as ratios of radii determined from the topol- 
ogy of the electronic density.^ Isoparametric transfor- 
mations then permits analytical mapping of polyhedra 
subdomains to a bi-unit cube, which are then simply in- 
tegrated by Gauss-Legendre quadratures. For any VP 
we only need evaluate numerically the integrands, Jaco- 
bian and weights at Gauss points for a relatively fast and 
accurate integral, and no issues with divergence. 

We showcase our isoparametric method by two means. 
First, we evaluate various electronic integrals analyti- 
cally using a well-known charge density model by van 
Morgan; 16 and show directly the accuracy and effi- 
cacy of our numerical method. Second, we implement 
the method in an all-electron Korringa-Kohn-Rostoker 
(KKR) code^ and apply it to a phase stability study 
of face-centered-cubic, (dis)ordered FePd. We exemplify 
the accuracy for formation enthalpies and the insensitiv- 
ity of results to the choice of spherical-harmonic basis-set 
(L-expansion) due to the use of weighted VP. 

We organize the paper as follows: After Section I 
background, we describe in Sec. II the RPC tessella- 
tion and weights that we merge with an isoparametric 
integration via an analytic dual-coordinate transforma- 
tions, known in the finite-element community, to cre- 
ate a general and optimal integration scheme. In Sec. 
Ill, we describe the van Morgan charge-density model to 
assess the performance of any integration method. In 
Sec. IV, we address the accuracy, convergence, and tim- 
ings of this isoparametric scheme for close-packed struc- 
tures; machine-precision accuracy with expected conver- 
gence is found, with millisecond timing for each VP. Our 
approach is greater than 10 5 faster and 10 7 more accurate 
than that with shape-functions. Finally, we discussed the 
results for application to FePd, then conclude in Sec. V. 
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I. BACKGROUND 

Given any atomic configuration, the first step to com- 
pute any site-centered integral quantity is to perform 
a Voronoi tessellation of the space in and around the 
molecule or solid, including possible empty sites to im- 
prove, for example, a site-centered basis set. Standard 
geometric (i.e., the Wigner-Seitz) tessellation subdivides 
space into VP such that every point within a VP cell has 
the property of being closest to one and the same site, and 
that each polyhedral face is orthogonal to and bisects the 
line segments joining the site centers, see Fig. [I] How- 
ever, in most materials (e.g., size-mismatched alloys), it is 
a poor subdivision of space, making VP corresponding to 
the smaller (larger) atoms too large (small) see Fig. [l] 
Such errors impact solid-state and biological problems, 
where, e.g., an accurate interstitial volume is required 
for reliable predictions of thermostability of cavity-filling 
mutants in proteins,^ or for statistical models in contin- 
uum systems where packing geometry plays a key roleP^ 

Given these deficiencies, various proposals have been 
made to place the dividing plane subject to atomic size. 
Richard® suggested using the ratio of the distance be- 
tween atoms and dividing plane to be equal to the ratio of 
the corresponding atomic radii, but this does not reflect 
bonding. RPC weights the distance to each atom by sub- 
tracting the squ ared atomic radius from the squared Eu- 
clidean distance j 13 * 14 ! which guarantees convex and space- 
filling VP, but radii must be provided. Other generaliza- 
tions include the introduction of non-planar boundaries 
between atomsPil For site-centered methods, the convex, 
space-filling property is critical; for example, in KKR the 
scattering matrices are only defined for convex VP, but 
the spherical-harmonic basis must reflect accurately the 
spherical density of each site or else the basis must be 
augmented (e.g., plane waves in FLAPWp^ or spherical 
waves^l). Notably, the VP tessellation benefits from a 
judicious choice of weights related to electron density.^ 

With a Voronoi tessellation in hand, a scheme must be 
chosen to perform VP integrations. In one-dimension, 
there are many numerical techniques that easily achieve 
accurate results. Yet, accurate and fast techniques 
for three-dimensional integrands with unusual domains, 
found in molecular or solid-state calculations, remain an 
active area of studyP^ We mention some key previous 
work. Ellis and Painted used Diophantine integration^! 
in molecular calculations, with convergence as 0(iV -1 / 2 ) 
to 0(N~ 1 ) for N sample points, just better than Monte 
Carlo. Becke^ used standard Voronoi partitioning with 
simplifications introduced to reduce multicenter integrals 
to a sum of single-center ones, with less detailed conver- 
gence studies than we provide here; we find that it uses 
slightly more points than ours for low (~10 -3 ) accuracy 
and much slower convergence for higher precision. More 
recently, Gaussian product formulas have been found use- 
ful when awkward domains of integration were split into 
tractable subdomainsPSlt is, however, not only the form 
of the integration domain but also behavior of the inte- 



grand that may necessitate the use of product formulas 
and further subdivision of the subdomains. 

For site-centered basis sets, there are limited choices 
of subdi visio n for convex VP domains. Baerends and 
coworkers^! broke each VP into subpolyhedra formed by 
the nucleus and its base of one of the VP faces followed by 
further subdivision of the face into connected sets of tri- 
angles and quadrilaterals. Averill and Painter 4 cropped 
each VP by an inscribed sphere to form an interstitial 
region associated with each VP face; hence, interstitial 
integrals are expressed as a sum of integrals over cropped 
pyramids. (In finite systems, a separate subdomain is 
the bounding part of the space outside the local atomic 
VP.) We use an analytic transformation to a bi-unit cube 
for the cropped integrals, an approach similar, but not 
the same as, Baerends and coworkers} 3 ^ as we discuss. 
Also, previous methods never considered using the un- 
derlying charge density to chose an optimal subdomain 
of integration. In all-electron methods, the wavefunc- 
tion and potential have cusps and singularities near the 
nuclei, respectively; these functions are easily integrated 
over spherical domains, although the "interstitial" region 
(between the sphere and VP facet) has a complex shape. 
Uniquely, we utilize the charge-density topology and the 
behavior of the integrands to determine the VP (spheri- 
cal and interstitial) subdomains. 

Although various integration methods have been pro- 
posed, the accuracy attained has usually been poor com- 
pared to the computational effort expended, typically 
a modest (7 — 8 digit) accuracy required a large num- 
ber of sampling points. The desire, of course, is to ap- 
proach machine-precision accuracy with a modest num- 
ber of sampling point s. W hile noting some similarities 
with previous methods)^ the present approach is unique 
in features and, particularly, its efficiency and accuracy 
for polyatomic systems; also it has the advantage of be- 
ing conceptually simple and easy to implement. We ver- 
ify that accurate (14 digit) integration over various ker- 
nels is achieved with a modest number of Gauss points. 
Moreover, when combined with the use of physics-based 
weighted VP, an insensitivity to site-centered basis set is 
possible while achieving high accuracy, as we show. 



II. METHOD 

The partitioning of space in and around a molecule or 
solid into convex polyhedra by RPC is described, followed 
by an analytic dual-coordinate transformation (isopara- 
metric mapping) of a bi-unit cube to obtain the shape 
of any specific facet subdomain of each VP. Dissection 
of each VP can be accomplished in two ways: Either 
each VP (1) is divided into subdomains formed by the 
pyramid between face and a sites origin; or (2) is split 
into an inscribed sphere domain and a sum of interstitial 
domains between sphere and each facet plane. We then 
find the integral at each VP as a sum over the Gauss 
quadratures by numerical evaluation of the integrands, 



FIG. 1: (Color online) Tessellation via (a) geometry and (b) 
RPC. The "size" of little (big) atoms is over (under) estimated 
in (a), affecting properties, whereas (b) reflects atomic "size" 
if radii are determined from charge densities, see text. 

weights and Jacobian of the transformation. If the inte- 
grand is evaluated at known Gauss points, only sampling 
points determine the error; whereas, if a numerical (dis- 
crete) grid is used, then interpolation error needs to be 
ameliorated. 



A. Weighted Voronoi Radical Plane Construction 

In unweighted tessellations, the subspace of points 
closer to an atom centered at than any other atom 
centered at r\,- is considered the VP about atom r^. For 
an atom on the edge of a molecule, the "polyhedron" 
will be unbounded; that can be amended by introducing 
an extra dividing plane that is tangent to the inscribed 
sphere of the unbounded cell. The best choice will be 
the plane that also minimizes the volume in the resulting 
cell. Of course, no such edge difficulties can occur for 
periodic crystals. In the case of a one-atom crystal, the 
VP is simply the Wigner-Seitz (geometric) cell. 

For a crystal with atoms of unlike size, the standard 
Voronoi partitioning of space assigns too much volume 
to the smaller atoms, see Fig. [TJa), and the correspond- 
ing VP is unphysicalP^ The biased partitioning incor- 
rectly determines the atomic volumes and the VP struc- 
ture. A weighted tessellation corrects by re-weighting 
the distances to the site centers. This re-sizing allows 
the volume given to each site to grow or shrink in accord 
with its size (radius). RPC uses the weighted metric 
||r — ri|| 2 — i? 2 , where Ri is the radius of the i th atom 
and ||r — is the Euclidean distance between point r 
and atom center rj. The advantage of the RPC tessella- 
tions is that the resulting domains are guaranteed to be 
convex polyhedral whose inscribed spheres match the 
input radii {Ri}, see Fig. [ljb). 

B. Physics-Based Definition of Atomic Size 

To provide the weights (radii) for RPC, we must choose 
the "size" (radii) of each atom. A simple choice is the 
atomic radii from empirical or theoretical tables, 27 which, 
however, is not site specific nor does it reflect bonding. 



FIG. 2: (Color online) VP for 2-atom unit cells with two 
inscribed radii. (Left) B2 cell with equal radii each of type A, 
and (Right) BCT cell with unequal radii of type A and B. 



For electronic-structure use, a judicious choice for each 
site (atomic or empty) are the minimum radii selected 
from the set of saddle-point radii (SPR) in the total 
electronic charge density, which reflects atomic "size" 
Initiating a calculation, these SPR can be estimated by 
overlapping the isolated-atom charge densities in the de- 
sired structural positions, similar to Lowdin potentials.^ 
For a spherical-harmonic basis, we have shown that this 
SPR representation, even within atomic-sphere approxi- 
mation, dramatically improves basis-set convergence and 
energetics in size- misma tched systems compared to full- 
potential methods! 15 ^ 28 ! Site-charges now also obey elec- 
tronegativity rules, as fou nd also with Bader topologi- 
cal (non-convex) cells E32H Full details of applications are 
available in Ref. IT51 

For RPC, given the site locations (structure) and SPR 
(weights), we use our modified version of Bernal's For- 
tran software^ to generate the VP information (vertices, 
faces, edges etc.). Figure [2] shows the VP generated for 
a B2 cell with atoms of equal size and a tetragonally- 
distorted BCT cell with atoms of unequal size. 



C. Dual Coordinate Transformation and Gauss 
Quadrature Sums 

Having divided the system into VP, there are two ways 
to proceed depending on the nature of the integrand f{r). 
For simple integrands, separate each VP integration over 
a numerous simple polyhedra associated with each VP 
face and perform Gauss quadrature sum, and the method 
works straightforwardly. If /(r) has singularities near 
the origin, or if it is accessible only on a sparse grid, then 
two major VP subdomains need to be handled separately, 
i.e. inside and outside of the inscribed sphere. If /(r) is 
spherical, the integral is one-dimensional and easy to per- 
form accurately, whereas the second, interstitial domain 
is more challenging, see Fig. [3] for FCC example. 

The interstitial has too unusual a boundary for the di- 
rect determination of suitable sampling points and their 
weights. To find the sampling points, we transform a 
bi-unit cube —1 < x,y,z < 1 into each pieces of the 
interstitial formed by each VP face and the site center 
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(a) 



FIG. 3: (Color online) (a) VP of a FCC structure with twelve 
quadrilateral faces and an inscribed (touching) sphere, (b) A 
section of the VP shown as single truncated pyramid. 





in Fig. |3jb), and introduce spherical coordinates (r, 6, 4>) 
so that the z-axis is perpendicular to the VP face. Within 
each piece, the radius r runs from the inscribed radius R 
to the pyramid base (or VP face). To consider the case 
where the inscribed sphere integral is not done separately, 
take R — > in what follows, and each pyramidal piece 
will no longer be cropped. 

Before we map the cube to this element, we must find 
a transformation that flattens the curved interior surface. 
Choose any three of the four corner vertices formed by 
the intersection of the pyramid and the inscribed sphere. 
These three points are taken to define an interior plane. 
Now consider a cross-section of the element at fixed angle 
<f> or 9, which resembles Fig.|4ja). Note l n is the distance 
from center of the inscribed sphere to point of intersection 
of radius vector with interior plane; and I / is the distance 
to intersection with base plane (or face). Then the map 



If — In 



l f (R-l n )+r'(l f -R) 



(1) 



(a) (b) 

FIG. 4: Cross-section of the cropped pyramid (a) before radial 
scaling and (b) after radial scaling. 

but cropped by the inscribed sphere. If any face has 
more than four vertices, points are added within the face 
(uniformly distributed) so that each face can be subdi- 
vided into polygons always having at most four vertices (a 
quadrilateral base); as a result, no interstitial subdomain 
has more than eight corners, like the cube. The same 
map used on the Gauss-Legendre points tells us the sam- 
pling positions in each interstitial subdomain. Note that 
one could use a triangular base, but we find that, while 
both subdivisions give the same results, the quadrilateral 
requires less operations, hence, it is more efficient. 

For clarity, consider a one-atom FCC crystal, as in 
Fig.[3ja), where the VP consists of 12 quadrilateral faces, 
which are divided into 12 cropped pyramids. Pick one, as 



will radially expand the interstitial piece (unprimed co- 
ordinates) so that the surface cut of the inscribed sphere 
will map to the interior plane (primed coordinates). Note 
that the map as given takes the plane to the sphere, be- 
cause, ultimately, we want a map from the cube to the 
interstitial piece. Despite the simplicity of the map (Eq. 
[I]), the Jacobian J\ is non-polynomial due to the angu- 
lar dependence of lf(6,<fr) and l n (9,4>). The standard 
determinant form of J\ can be simplified by considering 
the volume change of an infinitesimal cell embedded in 
a spherical coordinate mesh. The cell will be stretched 
radially by a factor of dr/dr' = (If — R)/(lf — l n ). And, 
because the cell will be translated radially from r' to r, 
the base area will change from r' 2 dfl to r 2 dfl. Thus, the 
total volume change (ratio) of the cell will be ^ ^ ■ 

Having flattened the interior, curved surface, we then 
perform a second mapping from this hexahedra to a bi- 
unit 2x2x2 cube, as depicted in Fig. [5j Let (a;', y', z') 
and (x" ,y" , z") be the coordinates before and after the 
transformation, respectively. Mathematically, we can 
connect them using the expression 
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where the index in the subscript (1 to 8) indicates the vertex number in Fig. [5] In this map, we have reverted to 
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FIG. 5: (Color online) Two-step coordinate transformation: 
(1) Bottom curved surface (a) to the interior plane (b) via the 
Jacobian Ji, and (2) hexahedra in (b) to the isoparametric 
(2x2x2) bi-unit cube (c) via the Jacobian J2. 



describe the hexahedral element in cartesian coordinates 
(x' , y' , z') rather than the spherical (r', 9' , </>'). 

The Jacobian of the transformation J2 that turns the 
hexahedra into a bi-unit cube is 



J 2 = 



dx dx dx' 

dx" dy" dz" 

dy' dy' dy' 

dx" dy" dz" 

dz 1 dz 1 dz' 

dx" dy" dz" 



(3) 



Thus, the volume integral over the interstitial region 
transforms to a volume integral over a cube. This can be 
expressed, using Gaussian-Legendre integration, as 

/ /(r)d 3 r = f C C dV'/(r") J x J 2 

JQIS J-lJ-lJ-l 



Ni N m N. 



l — l m — 1 n—1 

x wi(xi) w m {y'^) «;„«) (4) 

where J — J1J2, and Ni, N m and iV n are the num- 
ber of quadrature points along x"-,y"- and z"-axes, re- 
spectively. The Gauss points Xi and weights Wi are 
known analytically from the zeroes of the Legendre poly- 
nomial, so Eq. Q is straightforward to evaluate. Cal- 
culation time is primarily spent in numerically evalu- 
ating the analytically-derived Jacobians J\ and J2 for 
the two successive transformations and the f(x", y", z"), 
hence, quite fast. This isoparametric approach achieves 
machine-precision error for VP integrals involving vol- 
ume, charge-densities and potentials. The function 
f(x",y",z") should be evaluated at the specified Xi 
points; if, however, / is only defined on a discrete grid, 
the function must be interpolated to each Xi, in which 
case interpolation error is the major error that should be 
ameliorated to achieve high-accuracy integration. Gen- 
erally, if f(x",y",z") is a polynomial of order p\, p 2 
and P3 along the three directions, respectively, then 
the number of sampling points N required to integrate 
the quantity exactly for a simple polyhedra domain is 
+ 1) X (f + 1) X + 1). For the case where we sep- 
arate the integral over the inscribed sphere and integrate 
the interstitial over a domain that is curved, the transfor- 
mation makes the integrand effectively non-polynomial; 
therefore, more Gauss points will be required. 



Our method is distinguished from that in Refs. |3] and 
[5] by the choice of transformations, as well as partition- 
ing space via the weighted VP. We transform the Gauss- 
Legendre sampling points inside a bi-unit cube into the 
truncated pyramid by (1) cubic polynomial mapping of 
the corner points of the cube to the corner points of the 
truncated pyramid (given by J2), and then (2) perform- 
ing a linear mapping (in radius) of the interior plane (or 
side closest to origin) onto the relevant cut of the in- 
scribed sphere (given by J\). Our Jacobian J = J\Ji is 
always smooth and well-behaved, even for highly skewed 
pyramid. Baerends et alJ^ have noted that their choice 
of coordinates can cause their intermediate functions to 
behave poorly (i.e., the J diverges) when the pyramid has 
wide opening angles, or a strongly skewed face. For our J 
to diverge, the interior plane would need to (nearly) touch 
the pyramid base plane; but, with the interior plane de- 
fined as the one passing through three of the intersection 
points of the inscribed sphere and the edges of the pyra- 
mid, this could only happen if the sphere touched one 
of the corners of the VP, which can never happen. In 
addition, the present procedure requires minimally fewer 
function evaluations. 



D. Symmetry Considerations 

We could take advantage of the symmetry of the VP 
and crystal. We consider two kinds of symmetry. First, 
the point-group symmetry of the crystal structure iden- 
tifies the set of inequivalent sites in the cell, which re- 
duces computational time to that over inequivalent sites 
only. The second symmetry is associated with individual 
polyhedra. Because each polyhedra consists of various 
quadrilaterals associated with their faces, we can identify 
the symmetry equivalent pyramids by applying a set of 
symmetry operators over each polyhedra around the sym- 
metry unique atoms. By integrating only over symmetry- 
inequivalent pyramids corresponding to each symmetry- 
inequivalent atoms and weighting them with their degen- 
eracy, an appreciable savings in computer time would be 
obtained for systems with high symmetry. For example, 
there are 12 facets for VP in an elemental FCC structure, 
so, at a minimum, we could perform an integration over 
one VP facet and multiply result by 12; however, because 
each facet is 4-fold symmetric, we could do 1/4 of the VP 
facet and multiply result by 48, reducing VP integration 
by 1/48 of above integrations timings. 



E. All-Electron Implementation 

Our method is conceptually simple and easy to im- 
plement in any general electronic-structure code, with 
additional advantages for site-centered methods. For a 
spherical-harmonic basis, integrating separately over the 
inscribed sphere directly eliminates the Coulomb singu- 
larities due to the Jacobian within spherical coordinates. 
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Also, using the optimal SPR basis^we have better basis- 
set convergences and site charges. We have included this 
isoparametric method in a small set of RPC routines to 
determine the VP (vertices, faces, and edges). For com- 
plex cells, our algorithm has the flexibility to control the 
desired precision to balance the computational cost. 

This software is used to implement the (un)weighted 
VP-based isoparametric integration in our all-electron, 
KKR-CPA codeP^ We discuss these results in Section 
I1V1 Details of the calculations are as follows. The 
Green's function and related integrals use an external 
angular-momentum cutoff up to L max = 3, i.e., s— , 
p— , d—, and /—symmetries, as needed. Energy inte- 
grations of the Green's functions are done by contour 
integration^ via Gauss-Chebyshev methods with 18 en- 
ergy points. We use the local spin-density approxi- 
mation (LSDA) as parametrized by von-Barth-HedinF^ 
The Brillouin zone integrations use the Monkhorst and 
PacliP^l special k-point method with 20 3 points in the full 
zone. For disordered alloys, we use the coherent-potential 
approximation^ (CPA) based on the screened-CPA^ to 
incorporate more properly the metallic screening due to 
charge correlations in the local chemical environment. 



III. AN EXACTLY SOLVABLE MODEL 

To illustrate the numerical convergence and accuracy, 
we use Van-Morgan's exactly solvable charge-density 
modeli^ Many standard electronic-structure kernels can 
be exactly evaluated for the van Morgan density and po- 
tential, so the error in the numerical integrals can be 
precisely determined. We verify that accurate results are 
found with a modest number of Gauss points that depend 
on structure, and machine-precision can be achieved by 
increased number of points, slightly increasing computa- 
tional time. 

We showcase the convergence of volume and charge 
conservation, the [p(r)V(r)] integral evaluated for kinetic 
and/or Coulomb energy, and more highly varying func- 
tions in I and r. Apart from the cubic structures, we 
have also tested the convergence of the interstitial vol- 
ume integral for more complex crystal structures. In the 
timings below, we have not utilized the associated sym- 
metry of the crystal and the VP, so that the results reflect 
the most inequivalent case. 



A. The van Morgan Test Problem 

The van Morgan^ test charge density is defined as 

p(r) = i?fv T " r , (5) 

71=1 

where T„ are the nearest-neighbor reciprocal lattice vec- 
tors, and B is a scale factor. We will take B = 1 for 



simplicity. (From the Bauer expansion, a plane wave re- 
quires, in principle, an infinite number of spherical har- 
monics to be fully represented.) Because £l VP and il MT 
are known exactly for any crystal structure, it is often 
convenient, especially for site-centered methods, to di- 
vide the VP into two volumetric regions: the volume of 
inscribed sphere Q, MT and the volume within the inter- 
stitial region fl IS , so that Q vp = fl IS U Q MT . 

First, we can precisely assess the numerical error asso- 
ciated with volume conservation via 

/ d 3 r = n IS = n vp - n MT , (6) 

J is 

where fl MT = 4?ri? 3 /3, and, for example, R is 1/2, y/2/4, 
and V3/4 for SC, FCC, and BCC (in units of lattice 
constant), respectively. The left-hand-side numerical in- 
tegral is compared with the analytical result available for 
the right-hand side. For example, the VP volumes are 1, 
1/4, and 1/2 (in units of lattice constant cubed) for SC, 
FCC, and BCC, respectively. 

Second, we can assess the integrations associated with 
charge conservation, including the determination of elec- 
tronic chemical potential or Fermi energy. With p(r) 
having no zero-mode component in its Fourier expansion, 
the integral of charge over a VP cell must be identically 
zero; hence, charge neutrality requires that 

Q total = f p(r)d 3 r = 0. (7) 
Jvp 

Subdivision of VP yields 

Q IS = [ p{v)d 3 r = - [ p{r)d 3 r. (8) 

Next, we can assess numerical errors for the p(r)V(r) 
integral, which can be expressed as 

[ P V] IS = / p{r)V{v)d 3 r 

= - / P(r)V(r)d 3 r. (9) 

The exact analytic solutions of the above kernels for 
three cubic structures (SC, FCC and BCC) are given 
in the Appendix [X] However, non-cubic structures are 
numerically no more difficult or error prone than these 
cubic cases (but they cannot be performed analytically). 
Besides band-energy (an eigenvalue summation requir- 
ing a Fermi energy) and exchange-correlation, the above 
three integrals reflect the main integrations contributing 
to DFT total energies, for example. 

B. Complex Varying Integrands 

The general method we have presented here can inte- 
grate over any arbitrary polyhedra for any complicated 
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function, such as those with high angular momentum or 
strongly varying with(out) exponential decay. Here we 
showcase a set of strongly varying integrands that are 
critical for evaluating the near-field contributions to the 
Poisson equation (in the so-called moon region) for site- 
centered, electronic-structure methods,^ i.e., 



E 



drV(r') 
VP l r 



F im (r' + R) 



R| /+1 



(10) 



where Y\ m are the spherical harmonics and there is a 
Madelung summation over R. It is a rapidly decaying 
function with increasing and, to achieve high precision 
of this piece of the Coulomb potential, rather high I's are 
required. Convergence of the above integral for various 
{Zm}-values is shown in the next section. We have also 
tried other more strongly varying functions, and again 
achieved accurate results with modest number of Gauss 
points. As will be discussed elsewhere, most codes that 
implement the correct Poisson solution for space-filling 
VPs cannot do the integrals for skewed VP, or they are 
not accurate enough due to use of, e.g., shape- functions, 
an example appears below. 



IV. RESULTS AND DISCUSSION 
A. Accuracy 

To illustrate the convergence of isoparametric inte- 
gration, Figure [6] shows the logarithmic error in inter- 
stitial volume for six structures (i.e., 1-atom cubics, 2- 
atom hep, and 2-atom B2 and BCT). Each point on the 
graph represents the result for a combination of quadra- 
ture points (Ni, N m , N n ). From Fig. [HJa), it is clear that 
the cropped pyramid has a thinner dimension along the 
z-axis compared to the other two axes. Therefore, we 
use less quadrature points along z" than the x" and y", 
i.e., N n < (Ni,N m ); in particular, we used N t — N m . 
Accuracy of around 10 -3 is already reached with only 
Ni = N m = 4 points along the x" and y" . The darker 
line in each panel shows the minimum number of quadra- 
ture points along z" to achieve a convergence to 13 deci- 
mal places. For example, the minimum number of Gauss 
points along z" for a BCC structure to attain an error 
less than 10 -13 is two. The minimum number of points 
(N[,N m ,N n ) required is listed below each subpanel. 

The convergence of the charge density integral (Q) 
is given in Fig. [7j The left panel shows the logarith- 
mic error in the interstitial charge Q IS for the cubic 

structures. The right panel shows the absolute error 
c vp _ 



^calc 



QYJact m tne total charge integral. The 
charge convergence requires more points to yield a similar 
level of accuracy. For example, to achieve an accuracy of 
up to the third-decimal place, the BCC structure requires 
8-points along the i" and y" compared to the 4-points 
needed for the SC and FCC structures. Higher accuracy 




18x18x2=1176 




FIG. 6: (Color online) Logarithmic (base-10) error in inter- 
stitial volumes for six structures. Ni (i = I, m, n) is the num- 
ber of Gauss points along x" , y" and z" , respectively, with 
N„ < Ni = N m due to a smaller caliper along i" . Dark 
(blue) lines indicate minimum number of points along z" (to- 
tal points listed below plots) to achieve 13 decimal accuracy. 



requires more points for BCC case due to its wider and 
more asymmetric interstitial region. 

Figure [8] shows the convergence of the interstitial [pV]- 
integral for the cubic structures. It is interesting to notice 
that [pV] -integral converges almost at the same rate as 
the p-integral (left panel of Fig. [7]) and does not take 



TABLE I: Convergence for the interstitial volume, charge and 
[pV] integrals for various crystal structures. {Ni — N m ,N n } 
are the optimal number of points for each structure to reach 
an accuracy of at least 13-decimal places. VC, QC and [pV] 
stands for the volume, charge and [pV]-integral convergence. 

Structure {N h N n } vc {Ni,N n } QC {Ni,N n } pV 



SC 

BCC 

FCC 

HCP 

B2 

BCT 



{18,2} 
{15,2} 
{13,5} 
{12,5} 
{15,2} 
{12,5} 



{20,6} 
{26,8} 
{12,6} 



{18,8} 
{26, 10} 
{12,6} 




12x12x6=864 



FIG. 7: (Color online) For the van-Morgan problem for SC, 
BCC, and FCC, (left) the logarithmic (base-10) error in the 




18x18x8=2592 



BCC 




26x26x10=6760 



FCC 




interstitial charge, i.e., e = (Qcai, 



i exact 



)/Qexact> an d 



(right) absolute error in VP total charge, i.e., e = Q c 



' exact 



Other details are as in Fig. [fi] 



12x12x6=864 



FIG. 8: (Color online) For the van-Morgan problem for 
SC, BCC, and FCC, the logarithmic (base-10) error in the 
interstitial-[pl^] integral. Other details are as in Fig. [6] 



any longer for all the three structures. BCC structure 
in this case also requires comparatively more points to 
converge for the same reason in the previous paragraph. 
We have analyzed other relevant integrands (often used in 
electronic-structure calculation) as well and found either 
a similar or a slightly slower rate of convergence. 

In Tabic |TJ we have listed the minimum number of 
points required to get the interstitial volume, charge and 
[pF]-integral convergence to more than 13 th decimal for 
each structure. The number of points required are given 
as {Ni — N m , N n }. Table |ll] shows similar results for the 
full VP integral of strongly varying functions in Eqn. ( 10 ) 
for various {I, m} values, exhibiting oscillatory angular 
dependence with /-dependent spatial decay. In spite of 
its strongly varying nature, the number of Gauss points 
{TV} required to achieve an accuracy of up to 10-decimal 
place is not large, and are comparable to those of p and 
pV integrals shown in Table [TJ 

The accuracy of all our integrals is limited by the accu- 
racy of the VP boundary (vertices, faces and edges) infor- 
mation generated from the Bernal's softwarePS We have 
modified Bernal's original (binary-math/single-precision) 



code to improve its efficiency and extend its accuracy, and 
we were able to achieve just below 10 -13 . We have veri- 
fied that our main limitation in accuracy is due to lack of 
a double-precision real code. By rewriting the software 
from scratch, which is a considerable effort beyond the 
scope of present work, we could certainly achieve machine 
precision. Therefore, all integration results will be lim- 
ited to just below 10~ 13 ; with improved accuracy of VP 
information, machine-precision is achievable with similar 
Gauss points described. 



B. Efficiency 

To contrast the VP construction timings, we compare 
to the time required to expand the shape function (or 3- 
D step function) into spherical harmonics. 3 ^ The shape- 
function approach is often used in the community when 
needing site-dependent quantities. The EMTO, KKR, 
LSMS, APW, etc., methods, for example, typically re- 
ports site-quantities, and KKR Green's function methods 
require site-dependent VP scattering matrices. 
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TABLE II: Convergence of fee ai m in Eqn. ( 10 1 versus {I, m} 
(R is summed to 8th neighbor shell). {N} is the number of 
points per direction for 10-decimal place accuracy. 



/ 


in 


{N}a L 


\j^X Tfh ] num Brie til 


I Tn, ] exact 








12 


0.009951109455 


0.009951109341 


2 





12 


0.000000000000 


0.000000000000 


4 





14 


-9.449717387589 


-9.449717387292 


4 


4 


14 


-5.647286285886 


-5.647286285399 


6 





16 


-10.62648231314 


-10.62648231381 


6 


4 


16 


19.88032802119 


19.88032802174 


8 





1 Q 

18 


i?r o A C\f}A CI A A QC 


i?r o A C\C)A CI A A 1 n 

DO.o4U24ol441U 


8 


4 


18 


24.75927136401 


24.75927136434 


8 


8 


19 


37.72380770670 


37.72380770699 


10 





21 


135.8520785234 


135.8520785239 


10 


4 


21 


-136.8931058432 


-136.8931058438 


10 


8 


24 


-162.9353862901 


-162.9353862928 


12 





24 


-205.2050307885 


-205.2050307878 


12 


4 


26 


-493.6493920838 


-493.6493920829 


12 


8 


26 


544.8113627129 


544.8113627135 


12 


12 


26 


-261.8046415883 


-261.8046415883 



The shape-truncated function for a VP is denned as 

1 refi 



a(r) = 



r £ 



(11) 



where f2 is the VP region. The expansion of cr(r) in 
spherical harmonics yields the angular momentum de- 
composition 



a L (\r\) = £drY£(r) <r(r) = a L 



(r), 



(12) 



where the integration is over the angles f = (9, <fi) and 
L = (^m). The shape function is used to simplify the 
numerical integration of any function /(r) over the poly- 
hedron volume f2 as 

F = [ f(r)a(r)d 3 r 
Jn 

= fdrr 2 a L {r)f dtl Y L (r)f(r), (13) 

L=0 U 

especially if it is well-represented by spherical harmonics. 

The expansion coefficients <jl(x) must be truncated at 
a very high L trunc >> L max to achieve an accurate rep- 
resentation of the VP shape and to obtain a reliable in- 
tegral value. For example, for FCC structure, p(r) is 
well represented using L < 8 (i.e., L max = 8), but the 
shape- function should have L trunc » ^L max to have 
converged <jl<&{t) that will yield an accurate integral. 
As we shall see, this L trunc will limit the accuracy of 
the integrals in the codes that use this approach, mak- 
ing the shape-function approach unacceptable for general 
(non-high-symmetry) structures, where L trU nc should be 




10 12 14 16 18 20 



FIG. 9: (Color online) Timings to achieve a specific level 
of interstitial-charge accuracy for cubic structures using 
shape-function (left) versus isoparametric (right) integration. 
Shown in panels are logarithmic error in the interstitial charge 
(top), and times to construct VP boundary information (bot- 
tom) and to integrate (middle). Isoparametric integration is 
> 10 faster and achieves machine precision. 



significantly larger than in the cubic cases to achieve the 
same level of accuracy as FCC. 

Figure [9] shows accuracy and computer time for 
isoparametric (right panel) and shape-function (left 
panel) methods for SC, FCC, and BCC, for a direct com- 
parison. The rate of convergence is given with respect to 
the number of Gauss points along each dimension for the 
present method, and with respect to the l ma x for a fixed 
radial grid using shape-functions. The present method 
attains error in the van-Morgan interstitial charge below 
10~ 13 with less computational time. The shape-function 
technique cannot achieve an accuracy better than 10~ 7 
with l ma x = 16, an extremely expensive calculation due 
to the high-L expansion. Hence, our method provides 
some significant advantages over existing approaches. 

The bottom panel shows the time required to generate 
the boundary informations necessary to achieve a certain 
level of accuracy. For both methods, most of the time is 
spent in determining the VP boundaries. The present 
method generates this information in terms of neighbors, 
vertices, faces and edges for each VP. The shape-function 
method gets the VP shape in terms of an L-expansion on 
a specific radial grid. Clearly, the shape-function method 
requires > 10 4 more time than the present method. The 
middle panel shows the time (in msec) required to sum 
the final expression for the integration for both VP or 
shape function. The present method is faster by > 7 
times. Overall, using no symmetry (degeneracy) informa- 
tion to reduce the computational time, we achieve ~ 10 5 
faster integration with 10 6 less error. 
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TABLE III: Formation enthalpy AEf (in meV/atom) for 
(dis)ordered ferromagnetic FePd versus ~L max from equal (un- 
weighted) and SPR-weighted VPs, along with ordering ener- 
gies AE°- d . SPR results are much less sensitivity to basis-set 
L ranl cutoff. KKR results are compared to other results. 



TABLE IV: Excess charges within VP or ASA spheres for 
(dis)ordered FePd with (un)weighted VP via KKR. 



Lmax 


Unweighted VP 
AE° f rd AEf s AE°- d 


Weig 
AEf d 


;hted VP 
AEf s AE°~ d 


2 


-10.7 +18.6 -29.3 


-83.3 


-59.4 -23.9 


3 


-45.7 + 8.3 -54.0 


-88.8 


-63.9 -24.9 


4 


-30.7 - 9.9 -20.8 


-86.5 


-61.7 -24.8 




Expt. (Ref. EH 


-98 + 11 






VASP-LDA(GGA) 


+40(-130) 





C. Application to FePd 

We present the formation enthalpy, AEf, for FCC- 
based ferromagnetic (FM) (dis)ordered Fe-50%Pd, i.e., 
Al solid-solution and ordered Llo, from an unweighted 
(equal) VP and SPR-weighted VP, giving an optimal par- 
titioning of space for integration and concomitantly im- 
proved basis set, as described in Sec. |IIB| and shown in 
Fig. [2] Formation energies are highly relevant for phase 
stability; see Ref. EE for an example application to phase 
stability of magnetic-storage materials. 

Table [TTT] shows the AEf for FM Ll and Al phases 
le external L maa; for the local spherical basis. 



versus t 

Using the weighted- VP integration, results become sig- 
nificantly less sensitivity to L maa: , especially clear for 
the energy difference AE°~ d between ordered and disor- 
dered phases, which remain almost constant, in contrast 
to the unweighted case. Our weighted-integration yields 
formation energetics in very good agreement with that 
observed for Ll !^ The weighted- VP integration thus 
provides accurate results, a minimal basis set in terms of 
angular momentum cutoff, and a significant reduction in 
matrix-inversion time because of the now-permitted use 
of the lower rank of the KKR matrices N(L max + l) 2 , 
where N is the number of atoms in the unit cell. 

We could not find enthalpy data for the FM-A1 phase. 
Therefore, we provide AEp~^j in the paramagnetic phase 
(via disordered local moment state), which is related to 
order-disorder temperaturepSl Indeed, for the weighted- 
VP case, AE p -^ is 82 meV (or 952 K), close to the order- 
disorder temperature of 1050 Kj^ and showing that the 
disordered phase results are now accurate, too. 

Large-scale boundaries are critical in material science, 
e.g., for mechanical properties as pinning centers for 
mageto-anisotropy. As a test of the present weighted- VP 
integration, we calculated the [001] anti-phase boundary 
planar defect energy of Llo-FePd, a long-period structure 
with 32 atoms per cell with varying interstitial regions. 
We find 910 mJ/m 2 versus 890 mJ/m 2 from VASP plane- 
wave calculations, only ours is about 20 times faster, and 
provided local properties directly. 

In addition, the magnitude of site excess (or deficient) 
charge, i.e. the "charge transfer" , in a solid crucially 



Method 


Unweij 


dited 


Wei{ 


;hted 




AQ ord 


A Qdis 


AQ ord 


AQ dis 


VP Fe 


-0.111 


+0.059 


-0.032 


-0.026 


Pd 


+0.111 


-0.059 


+0.032 


+0.026 


ASA Fe 


-0.139 


+0.089 


-0.082 


-0.051 


Pd 


+0.139 


-0.089 


+0.082 


+0.051 



depends on the way in which the space is divided into 
geometric cells. For space-filling VP, SPR-weighted cells 
will provide a more physics-based partitioning of space 
and more realistic assessment of charge transfer. Ap- 
proximate methods like the popular atomic-sphere ap- 
proximation (ASA) has an overlap error; the situation 
becomes worse for non-close-packed materials. 
In Table llV 



we show the calculated excess charges 
within (un)wcighted VP sites in Al and Ll FePd, with 
comparison to (un)weighted ASA spheres used in many 
popular codes. Reference [T5] provides details of ASA ap- 
proach. Generally, there is a charge transfer from Fe to 
Pd, as expected from the electronegativities. 

However, for unweighted cases in the Al phase, there 
is an excess charge on small (Fe) atom, distinctly un- 
physical, and due to the tails of the charge density of 
large (Pd) atoms being improperly cut off at the smaller 
radii. When a weighted VP or ASA is used, this situation 
is corrected (the sign changes) because the charge den- 
sity is now better represented in the disordered phased 
For the unweighted Ll case, there is a large depletion of 
charge on small (Fe) atom due to a Madelung effect; how- 
ever, for the weighted case, the inscribed sphere reflects 
more appropriately the extent of the charge density and, 
hence, it is a more reliable estimate. Importantly, there 
is a large reduction in excess (depleted) charges for the 
weighted- VP integration compared to the weighted-ASA 
case (now with the correct sign), which shows the error 
arising from overlap of spheres. 



V. SUMMARY 

We have presented a fast, accurate, and easy to imple- 
ment method for the numerical integration over general 
VP for polyatomic systems. The algorithm combines a 
weighted Voronoi partitioning of space with isoparamet- 
ric integration using the Gauss- Legendre quadrature for- 
mulas of product type, and does not suffer from any ill 
behavior with shape of VP. In contrast to other meth- 
ods, accuracy and convergence was tested rigorously via 
an analytic charge-density model, with machine-precision 
accuracy for reasonable number of Gauss points. We 
showed also that our algorithm is 10 5 faster and 10 7 
more accurate than that based on shape-functions used 
in several electronic-structure codes. Our method could 
be used for other types of condensed matter problems 
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requiring integration over arbitrary convex VP. Here, 
we implemented the general method in an site-centered, 
electronic-structure code and calculated formation en- 
thalpies for FePd, yielding good agreement with experi- 
ment. The radii to set the Voronoi/Delauney tessellation 
weights is obtained from a physics-based definition, i.e., 
the saddle-points in the total electron density. 

Acknowledgements: Our work was supported by 
the U.S. Department of Energy BES/Materials Science 
and Engineering Division from contracts with Illinois 
(DEFG02-03ER46026), Ames Laboratory (DE-AC02- 
07CH11358), operated by Iowa State University, and 
Lawrence Livermore National Laboratory (subcontract 
B573247). DDJ and SNK acknowledge support from 
NSF (DMR-07-05089) and, recently the Center for De- 
fect Physics, Energy Frontier Research Center. We thank 
Wan-Tsang Wang (NSYSU, Taiwan, and advisor Y.C. 
Chang) for help in algorithm testing while in Illinois. 



where a<» = TW R { a % , T< s ) = |T„|. is a normal- 

ization constant whose value is 24, 32 and 48 for the SC, 
FCC, and BCC lattice, respectively. Finally, the exact 
expression for pV integral ^ for the van Morgan charge 
and potential for any general lattice can be simplified for 
the cubic lattices (SC, FCC, and BCC) as 



exact 



[ T f.s)]2 



(A3) 



4?r 



[T( s )] 3 



(s) (s) (s) 

sin 7> — 7> cos 7. 



[A 



where = K = 6, 8, 12 and = 2, 3, 4 for the three 
lattices, respectively, and 7^ = /3j a^ s \ 

For the above integral expressions, the coefficients fi 
and Pi for the SC, FCC, and BCC lattices are 



Appendix A: Exact Solution for Cubic structures 

For volume conservation the VP and inscribed volume 
are known analytically, the interstitial integral is 



IS 



fi 



VP 



in 



t R (s) 13 
1 MT' ' 



(Al) 



for SC, A = -p— =K ; A = s/2/3 2 = 2. 



for FCC, h 



K 2 2 

K-2- 



K -2 



h = K; 



(3i = V2p 2 = \/-P 3 = 2. 



where (s) indicates the lattice type (SC, FCC or BCC). 



Rjmt * s t ne mscr ibed sphere or (muffin-tin) radius of the 
lattice s. For charge conservation, it is straightforward 
to show that 



^ exact 



A^TT 
[TW]3 



(A2) 



for BCC, fx = ^- A h = j^h = \h = K- 



Pi 



V3 



(3 2 = 2/3 3 = V2/3 4 = 2. 
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